DYNAMICS NEAR A MINIMAL-MASS SOLITON FOR 
A KORTEWEG-DE VRIES EQUATION 



J.L. MARZUOLA, S. RAYNOR, AND G. SIMPSON 

Abstract. We study soliton solutions to a generalized Korteweg 
- de Vries (KdV) equation with a saturated nonlinearity, following 
the line of inquiry of the authors in [5] for the nonlinear Schrodinger 
equation (NLS). KdV with such a nonlinearity is known to possess 
a minimal-mass soliton. We consider a small perturbation of a 
minimal-mass soliton and identify a system of ODEs, which mod- 
els the behavior of the perturbation for short times. This connects 
nicely to a work of Comech, Cuccagna & Pelinovsky [l]. These 
ODEs form a simple dynamical system with a single unstable hy- 
perbolic fixed point with two possible dynamical outcomes. A par- 
ticular feature of the dynamics are that they are non-oscillatory. 
This distinguishes the KdV problem from the analogous NLS one. 



1. Introduction 

We consider a generalized Korteweg-de Vries equation equation of the 
form 

(1.1) u t + d x (f(u))+u xxx = 

where / is a saturated nonlinearity; that is, / behaves subcritically at 
high intensities and supercritically at low intensities. An example of 
such a nonlinearity is 

(1.2) /(*) 



1 + 5sP~i 

with p > 5 and 1 < q < 5, and with 8 > as an additional parameter. 
In the computations we present, we take p = 6, q = 3 and 5 = |. 
A traveling wave function u(x,t) = c (x — ct), is a soliton solution 



to (1.1) when the profile <p c satisfies the ODE 
(1.3) - c<p c + f{<f> c ) + d yy 4> c = 0. 
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Figure 1. M c = H^dlia possesses a local minimum for 
the saturated nonlinearity given in (1.2). 



If we were considering a noncritical power nonlinearity, f(s) = s p , the 
equation would admit solitons of arbitrarily small L 2 -norm. However, 
saturated nonlinearities are more restrictive. For the nonlinearities we 
consider, there will be a unique soliton of minimal L 2 -norm. See Figure 



[T] for a plot of the L 2 -norm as a function of c for one instance of (1.2). 
This critical soliton is denoted (p Cir . 

It is known that there exists an interval U C R so that there is 



a soliton solution to (1.3) for each c G U |3j; moreover, 4> c (y) is a 
smooth function of c on U. Indeed, via elliptic theory we see that it 
holds generically that <fi G C q+2 ; since <fi > 0, this then implies that 
<f) G C°° for p, q G Z. Solitons can be interpreted as minimizers of the 
Hamiltonian energy 

(1.4) E(u) = \f \d x u\ 2 - I F(\u\)dx, 

where F(s) = Jq f(t)dt, subject to the fixed momentum condition 

1 r. l2 . 



1.5) N(u) = - I \u\ 2 dx, 



with c acting as the Lagrange multiplier. Equation (1.1) also conserves 
the mass: 

(1.6) I{u) = / udx. 

Jr 

By evaluating these conserved quantities at the soliton C , we get func- 
tions of c: 

(1.7) S c = E(<f> c ), K = N(<f> c ), X C = /(0 C ). 

The minimal-momentum soliton is found at the value of c = c* such 
that 

(1.8) ^K= (0 C ,<9 C C > =0. 

dc 

This first order condition would also hold at a maximal momentum 
soliton. 



While the stability of soliton solutions of (1.1) is well understood 
away from such critical points, our understanding of the dynamics near 
this critical soliton remain incomplete. Due to the saturated nature of 
the nonlinearity, the equation is globally well posed in H 2 , so there is 
no finite-time singularity. But does the perturbed soliton converge to 
some nearby stable state, disperse, or engage in some other dynamic? 
It is known [I], that the minimal-mass soliton itself enjoys a purely 
nonlinear instability. The purpose of this work is to further examine 
the dynamics of this type of solution. 

To better understand these dynamics, we consider perturbations of 
4> Cic . Beginning with the ansatz 

(1.9) u(x, t) = 4> c (x — ct) + p(x — ct, t) 

for a perturbed soliton, we obtain the following evolution equation for 
p: 

(1.10) p t = d y [-d 2 y p + cp- f{<p c )p] + 0(p 2 ). 

In order to analyze this equation, we first consider, in Section [2j the 
spectrum of the linearized operator, A c , where 

(1.11) A c = d y L c , L c =-d yy + c- f(<t> c ). 

An examination of A c reveals that its generalized kernel has a dimen- 
sion of at least two. At c*, this dimension increases to least three, and 
could be four under special circumstances. Thus, there will be secular 
growth of the perturbation at a critical point generated by the compo- 
nents of the solution parallel to elements of the generalized kernel. Even 
if the perturbation is initially orthogonal to these unstable directions, 
the higher order terms will likely generate unstable contributions. 
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Such secular growth is eliminated by making a more general ansatz 
that allows c to modulate in time. This is given in Section [3j where we 
introduce a three dimensional set of scalar parameters, including the 
soliton speed and wave center, and allow them to modulate. This per- 
mits for projection away from the linearly unstable modes. We separate 
the projection of p onto the discrete spectrum of A c from its projection 
onto the continuous spectrum. We then discard the continuous projec- 
tion component while modulating the remaining parameters to obtain 
a two dimensional system of first-order ordinary differential equations. 

Finally, in Section 111 we describe the numerical methods we use to 



compute key parameters and simulate (1.1). Section ^ presents the 



results of these computations and discusses their implications. 

2. The Linear Operator 
In this section, we survey the spectral properties of A c as defined in 



(1.11), about c = c*. For all values of c that admit a soliton, one can 
directly compute 

(2.1) A c {-d y (j> c ) = 0, 

(2.2) A c d c <p c = -d y <j) c . 

For the remainder of this article, we will denote differentiation with 
respect to c by '. The first two elements of the generalized kernel of 
A c , then, are: 

(2.3) ei )C = -dy<j> c , e 2 , c = <\> c - 

At a minimal-momentum soliton (or indeed, any soliton satisfying the 



first order condition (1.8)), there is a third independent function e^, 



in the generalized kernel of A c , which satisfies 
(2.4) A^e 3 ^ = e 2 „ 



To see why such a state exists, consider the adjoint operator, A* 
—L c d y . We immediately compute 

(2.5) A* c <p c = 0, A* c D- x <f/ e = -0 C , 
where 

(2.6) D- l f = [ V f. 



Hence we define: 
(2.7) g x , c = C , g 2 , c = D 
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Then, by the Fredholm alternative, we see that A c f = e 2iC has a solution 
provided 



(2-8) e 2 , c9l , c = I U = ^r f Ul = K 



dc / 2 



vanishes, which is condition (1.8). Thus, at a minimal-mass soliton, 
there is indeed a third element of the generalized kernel of A Cii , which 
solves 

A Ci ,e 3 ^ = e 2 , c ,. 

Consequently, there is also a third element in the generalized kernel of 
(2-9) Alg 3 ^ = g 2 ^ 

To see if there is a fourth element in the generalized kernel, we again 
consider A^f = e 3jC>t and recognize that we would need 

= (e 2 ,c,,^2, c J = J (p' c , (J 4>'cA dy 



(2.10) 



2 

I 

dy2\J_ OQ ~ rc *J 2 vX - 



d 1 / r y \ 1 



to vanish. Even at a critical value of c* for which Ml = 0, it is not 
generic to observe X' c = 0. Indeed, for the particular nonlinearity / 
that we consider, our minimal-momentum soliton will not have this 
fourth element. It would be of interest to find a nonlinearity that does 
satisfy this additional degeneracy condition, and to study the dynamics 
near the resulting doubly-critical soliton. 

A particular challenge, discussed below in Section [4], is that some 
of these generalized kernel elements, notably e 3)Cjt , # 2 ,c* an d <?3, Cy , are 
not in L 2 . While vanishes exponentially fast at +oo, it is only 
bounded at — oo. (?2,c* and both vanish at — oo, but at +oo the 
former is only bounded and the latter grows linearly. The reader can 
find a discussion the function spaces in which these kernel functions lie 

CD- 

For later use, we remark that away from c*, using the implicit func- 
tion theorem, there is a scalar A c and function e3 jC , both smooth in c, 
such that 

(2.11) (A c -A c /)e 3 , c = e 2 , c , \ c = - M * 

(0c,e 3 , c ) 



3. Modulation Equations 



To overcome the secular growth due to the generalized kernel, we 
now permit the equation parameters to modulate about the extremal 
soliton. First, we define the moving frame 

(3.1) y(x,t) = x — / c(a)da — £(£). 

Jo 

Next, we consider a solution u(x,t) which is a perturbed, modulating 
soliton, 



u(x, t) = (j) c (t) [ x - J c(a)da - £(t) 



(3.2) 

+ p[x— I c(a)da — £(t),t 



My) + p(y,t). 



Substituting this into (1.1), we get 

(3.3) p t + (f)' c c - d y (f) c i - p y i 

= d y L c{t) p - ^d y {f"(<p c )p 2 ) + higher order terms. 

We let 

(3.4) rj(t) = c(t) - c„ 
and decompose 

(3.5) p(x, t) = C(t)e 3 , c + v = C(*)e 3)C . + C(t)v(t)e' 3 ^ +v + 0(( V 2 ) 
as in [l], equations (3.19-3.24), which results in 

v t + A c v = -£e 1>c - (?) - C)e 2 , c - (C - ^cC)e 3 , c - f]d x p + d x N. 
To close the system, we introduce the constraints 

(9i,c,v) = (g2, c ,v) = (g 3 ,c,v) = 0. 

We make several observations about the resulting dynamical system. 
The term d x N has quadratic and higher-order terms. We will preserve 
only quadratic terms in our computations. We will also disregard all 
coupling to the continuous spectrum, though some of this may also be 
of quadratic order. Thus the quadratic order terms we are considering 
in our approximation can be taken to be ^( 2 d x (f" (4> c )e^ c ) . Projecting 
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onto the canonical spectral functions, the finite dimensional system 
then takes the form 

/ i \ ({9i,cA, c -\Uf"{<i>c)elc)Y 
(3-6) i) - C = -CX(C)- 1 (92, c , e 3 , c - \d x {f" (0c)ei, c ) > 

Under this approximation, £ is slaved to 77 and (. There is only weak 
coupling between £ and the other parameters through v. 
Continuing with the above assumptions, 

(— (<7i,c> 9 x e 3:C ) (<?i, c , e 3c ) 0\ 
-{92, c ,d x e 3 , c ) (g2, c ,eQ = % + (S c , 
-(93,c,d x e 3:C ) <^3, c ,e 3)C ) 0/ 

where (7^)^ = (gj, c ,^k,c)- Assuming 7^ and S c are 0(1) and ( is suffi- 
ciently small, we can make the approximation <S C ~ %, and thus obtain, 
for appropriate vectors R c and Q c : 

f(9i,c,e' 3 , c -ld x (f(<j> c )el c )y 
= -eS^y 1 <<? 2 , c ,e 3)C - \d x {f"{<t> c )el c )) 
\(fe,4, c -^,(/%)ey} y 

/<^c,e 3 , c -|^(r(0 c )ey>\ 
= -CX' 1 {92,cA,c~ \d x {f'(<t>c)el c )) +0(C 3 ) 

= -c 2 t c - 1 j r c + o(c 3 ) 
= -c 2 g c + o(c 3 ) 

Therefore, the leading order equations, subject to these approxima- 
tions, are 

(3.8a) i = -Q ljC C 2 , 

(3.8b) 77 " C = -Q 2)C C 2 , 

(3.8c) C - A C C = -Q 3 ,cC 2 - 

Making the Taylor expansions of A c the Q? )C 's about c = c*, and omit- 
ting the £ equation, we obtain the quadratically nonlinear ODE system 
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(3.9a) fj - C = -Q 2 , C *C 

(3.9b) c - KM = ~Qs,cX 2 - 

Critical points of the system are found at ( = 0, r\ = i] where i]q is 

arbitrary, and at £ = r] — x ,q 2 ■ The latter isolated critical point 

is a saddle point in the first quadrant. 
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Figure 2. A plot of a representative phase plane dia- 
gram. The hyperbolic fixed point is indicated by the o 
while the nonisolated critical points are indicated by the 
solid line ( = 0. 



For the ( = critical points the linearized problem is 

1 

The eigenvalues are A c+ ?7o and 0. Thus, depending on the sign of A' Ci ?7o, 
the solution is either linearly stable or unstable. It is worth noting 
in this context that A c = ,7^ c \ . Thus, A' = 7^ c * , since Ml = 0. 
Also note that, according to fl], near c* we have that (0 C , e3 jC ) > 0. 
Thus, we expect that the sign of A' c+ depends solely on whether we 
are at a minimal or maximal soliton. In each case, we expect to see 
that, depending on the sign of T]q, the critical point at (0, r] ) is either 
linearly stable or linearly unstable. There is semi-stability at (0,0), 
depending on the sign of the initial perturbation rjQ. See Figure [2] for 
a representative saturated KdV phase plane diagram. 



4. Computational Methods 



In this section we briefly outline the computations needed to make 
a comparison between (3.9) and the KdV equation, (1.1). Motivated 
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by (3.2), (3.4) and (3.5), we will take initial conditions of the form 



(4.1) Mo = 0c* (x) + r/ e 2 , c , + Coe 3 , c ,- 

4.1. Spectral Computations. To compute the coefficients appearing 



in (3.9), the Q/s and A^, we must compute 

• The generalized kernels of A Cir and A*^ from which we can get 
the inner products (<7j,c*> e fe,c*) which the matrix T Cic comprises; 

• d x e^^ and e' 3 ^ to obtain the R Cir vector, which, with , allows 
us to obtain the Q Ct vector; 



• A' c , which can be obtained by differentiating (2.11) and then 
computing <f)". 

The first few elements of the generalized kernel, e 1)Cii = d x cj) Ct , e 2)Cyi = 
(f)' and gx >Ci! = (f) Cir are readily obtained using the the sine discretization 
method previously used by the authors in |'5||. Briefly, this approach 
solves equations like 

(4.2) LJ = g 

using a sine discretization of /, g and L c , provided g is localized. 
Derivatives of functions are easy to obtain using the discretized dif- 
ferentiation matrices, and L 2 inner products are just finite dimensional 
inner products, multiplied by the uniform grid spacing. To obtain the 
minimal-momentum soliton, we use a root-finding algorithm to solve 
M' c = 0. See j4||6||8j[9] for additional details on the sine discretization 
method. 

Computing the other elements requires a bit more care as they are 
not L 2 -localized. First, let 

(4-3) Q(y) = I" e 2 , Ci . 

Then e 3jCit solves 

(4-4) L c ^e 3 ^ =6, lim e 3 ^(y) = 0. 

Given sine discretization of <fi' c , we can readily integrate to obtain 
using the techniques given in |tj. The use of such a quadrature tool 
was not required in [5] as all of the NLS kernel elements belonged to 
L 2 . This also gives us <?2,c*- 

Assuming that e^ tCji grows, at most, algebraically at — oo we can drop 
the /'(^cj^c* at large negative values of y term to estimate 

(4-5) - d 2 x e 3 ^ + c,e 3 , c . « -2^ + 0. 
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This approximation implies that is actually bounded at — oo: 

(4.6) lim e 3 , c , = -i^. 

j/-»~oo C* 

One method of computing e^ Cir is to split it into a piece which has the 
above asymptotics, and a spatially localized piece, 

(4.7) es = e ( 1 ) + e W 

where e 3 ^ G L 2 and ef^ vanishes at +00. Using the above estimate, 
we set 

(4.8) 1 1 

= --2^.-(l+tanh(-x)). 

We have some flexibility in selecting x~- The essential feature is that 
it should not contribute anything at +00, while capturing the known 
asymptotic behavior at —00. We then solve 

(4-9) L c ^l = Q-L c J*l 

As the righthand side is now localized at both ±00, we obtain e^. 

The adjoint problem is similar, but requires slightly more care. First, 
we solve 

(4.10) ^c>3,c, = <?2,c* 

and then integrate to obtain . While e^, Ci , wa s asypmtotically con- 
stant, g^ Ci< will have linear growth at +00. The other function which 
requires such an asymptotic splitting is e' 3 Cir , 
To compute A' c , we compute 



(4.H) A' = 



<f>'^ is /^-localized and obtained by solving 

Computing the various inner products, we obtain the matrix T c „ and 
the vector R Ct , from which we can solve for Q Ci< . This provides us with 
all coefficients in the ODE system. 
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4.2. A Finite Difference Method for KdV Type Equations. In- 



tegrating (1.1) with initial conditions of the form 
( 4 - 12 ) 0c*+m, + Coe 3)Cyt 

requires some care, as e^ tCir is not spatially localized. However, since it 
is asymptotically constant, we can, to leading order, use approximate 
Neumann boundary conditions at ±a; max , the edges of our computa- 
tional domain: 

(4.13) d x u = d 2 u = 0, at x = ±x max . 



Using this approximation, we can then solve (1.1) using a linearized 
implicit method formulated in [2]. This method has second order ac- 
curacy in time, with spatial accuracy given by the quality of our finite 
difference approximations of the derivatives. In this work, we use sec- 
ond order symmetric estimates of the first and third derivatives. 

4.3. Extrapolating and Matching Discretizations. A challenge 
in using our sine approximations of the kernel functions is that they 
are given on one discretized mesh which may not be sufficiently large 



to employ the approximate boundary conditions (4.13) for our time de- 
pendent simulation. To overcome this, we use the farfield asymptotics 
of these elements to extrapolate onto a larger domain with a given mesh 
spacing. 

Numerically, we discretize on a short interval, [— -R so i, -Rsoi] to com- 
pute the soliton using the iterative sine method from [HJ. We then 
asymptotically extend uq to a much larger interval, [— _R as ] for R as 
large relative to where we desire to have the boundary. In particular, 
we extend using the asymptotics 

0c* 0*0 ~ a\t~^ x \ as x ->■ ±oo, 

e 2,cA x ) ~ a 2 xe~ v/5 *' :r '', as x —7- ±oo, 

e 3,c*( x ) ~ a 3 x 2 e~^^, as x — > oo, 



e3,cA x ) 



-d c T - a A x 2 e~^ x \ as x ->■ -oo. 



The asymptotics of are standard, those of e2, c „ arise from the commu- 
tator relation [A,x]f = 2d x f, and those for e^ tCi , arise from integration 
by parts, given the asymptotics of e2, c *- To observe that such a continu- 
ation is nicely continuous and avoid boundary effects from the iterative 
methods, we actually choose to extend 0, e2, c * and from values de- 
termined of distance 1 from ±R so i. We include a log plot of an initial 
condition with i] Q , ( > in Figure [3} We then a linear interpolation to 

have an evenly spaced grid on another still large interval, [— R p d e , Rpdc] 

11 



10' 



10" 



10" 



s 10 



10" 



10 



10 



-100 



-50 



50 



100 



Figure 3. A log plot of the extended initial data. 



with Rpde < Rus, on which the simulation is perofrmed. Here, _R p d e 
is chosen large enough to minimize boundary interactions during the 
numerical integration, which in KdV simulations appear quite quickly 
due to the dispersion relation; see Figure |4j 

To summarize we have selected three different domains, i? so i , -R as 
and -Rpde, on which we respectively compute the soliton, match asymp- 
totics at ±oo, and then linearly interpolate onto a uniform grid to 
solve the PDE, with R so \ < R p d e < -R a s-The simulation then proceeds 
with a linearly implicit finite difference scheme using a split step time 
discretization. In our results, our PDE simulation contains many os- 
cillations that diminish as the boundary effects are minimized. 

4.4. Extracting Parameters. As we aim to compare the ODE sys- 



tem (3.9) with the PDE, we will need to find a way to extract £, 
r] = c — c*, and ( from u(x, t), 

(4.14) u(x, t) = c (x - / c - + Ce 3)C (x - f c- {) +v(x - f c- {). 

We estimate the wave speed by computing the center of mass of u at 
each time step, and estimating its speed by finite differences; this gives 
us c = c* + rj. Unfortunately, there is some ambiguity between the rate 
at which the wave moves due to the speed, c, and changes in the phase, 
£; we are only able to estimate, collectively, 



(4.15) 



d t y(x,t) = - c -£, 
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Figure 4. A semilog plot of the numerical solution to 
our PDE at t — 1 showing the oscillatory back scattering 
from the boundary at x = R. 

and assume that this is dominated by the c, at least for the time scales 
we study Indeed, on the time scales over which we numerically inte- 
grate, the £ term is quadratic in (, which must remain small for our 
computations to remain accurate. We then integrate the wave speed 
by quadrature to estimate the shift. 

Next, we estimate ( by projecting onto g x \ 

(u(-+,t),gi, e ) - \c{t) 2 {eL c ,g liC ) 



(4.i6) at) 



(e3,c,9i,c) + c*(e' 3 ,,5-i, c ) 



which is done by assuming that on the time scales we consider g l c , e' 2 , 
e3 jC and e' 3c are well- approximated by their values at c = c*. 

5. Shadowing Results 



We now take as an initial condition (4.1) and study the evolution for 
different r] Q and Co- 

The results appear as Figures [5] and [6j in which we take initial data 
that begins in the first and third quadrant of the phase plane [2] and 
compare the projection of our integrated numerical PDE to the pre- 
dicted ODE dynamics with domain size R p d e — 120.0, iV = 10 6 spatial 
grid points, time of integration T = 30.0, and time step h t = 10~ 4 . The 
remaining cases display rather similar behaviors. In [7], we observe that 
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Figure 5. A comparison between the PDE and the 
ODE with when r] = 5 x 10~ 4 and Co = 5 x 10~ 4 . 





Figure 6. A comparison between the PDE and the 
ODE with when r] = 5 x 10" 4 and ( = -5 x 10" 4 . 



by taking larger initial t]o > 0, Co < 0, our solutions diverge from the 
predictive dynamics on a shorter time scale (T = 20.0) with otherwise 
comparable parameters as above. Since the nature of KdV is to move 
to the right, in order to lessen boundary interaction, we solve the PDE 
in a moving reference frame around base velocity c* and implement the 



schemes to project onto the ODE parameters c and £ as in Section 4.4 



Implementing these approximative schemes and comparing the evo- 
lution of the corresponding ODE in (3.9), the figures show that for 
long times the dynamics indeed fit the predicted dynamics. For small 
enough perturbations of the minimal-mass soliton, the dynamical sys- 
tem predicts that the orbits travel very slowly towards the stable or 
unstable manifolds. Hence, we only follow the orbits on time scales 
where the parameter ( has made a large motion in its orbit. The c 
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Figure 7. A comparison between the PDE and the 
ODE with 770 = and Co = -10~ 3 . 

component varies essentially linearly on this scale however. As a re- 



sult, we observe that (3.9 ) is a good model for perturbations close to the 



minimal mass. As in the NLS case, there is a perturbation of the min- 
imal mass solution leading to dynamics that move c to smaller values 
where the corresponding solitons are linearly unstable. However, con- 
tinuing along this trajectory is forbidden by mass conservation as seen 
in Figure [1} We postulate, as we did for the corresponding Schrodinger 
dynamics in |5|, that this could be an energy transfer mechanism to 
the continuous spectrum in the infinite dimensional system. This would 
eventually lead to dispersion. However, as we are here working with 
perturbations that are not in L 2 , it is not possible to compare to known 
dispersive solutions as was done for Schrodinger. 

Appendix A. Details of Numerical Methods 



Using the sine methods described in Section 4.1 and similarly applied 



in 5 6, we compute the parameters for system (3.9). The convergence 



of these parameters, as a function of the number of grid points, is given 
111 the Table ft 

The kernel functions are given in Figure [8j As discused, e^ Cir , #2,0* 
and g3 jCi are clearly not in L 2 . 
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Table 1. Convergence of the ODE system parameters 
as a function of the number of grid points, M. 



M 


c* 




a; 


Qi 


Q 2 


Qs 


20 


0.76419938 


-0.18921573 


4.90659351 


4.14952911 


-1.21924052 


40 


0.76214845 


-0.19352933 


3.33351916 


5.02642785 


-1.35439215 


60 


0.7621^ 


^663 


-0.19276719 


2.51228473 


5.22964908 


-1.36903497 


80 


0.7621J 


J815 


-0.19262966 


2.21019564 


5.28041633 


-1.37169953 


100 


0.7621^ 


^822 


-0.19260110 


2.10124978 


5.29465657 


-1.37230370 


200 


0.7621^ 


^823 


-0.19259143 


2.03532641 


5.30133115 


-1.37253316 


300 


0.7621^ 


^823 


-0.19259139 


2.03442654 


5.30138737 


-1.37253443 


400 


0.7621£ 


^823 


-0.19259139 


2.03440251 


5.30138854 


-1.37253445 


500 


0.7621^ 


^823 


-0.19259139 


2.03440202 


5.30138858 


-1.37253445 
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Figure 8. The elements of the generalized kernels of 
A Ct and A Cii *, computed with M = 500 grid points. 
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